function [ssolve_market_all, ssolve_market, nnon_converge_all, nnon_converge]...
    = cf_rate_adjustmet_cost_lapse(lapse_u, medi_cf, scale_vec, cost2_0, cost2_1_all, D_input, D, N,...
    D_major, arg_share, D_nF, N_st,...
    D_enroll_st, ...
    alpha, ...
    Delta_xi_all, cost1_e_all, mu_f_all, sigma_f_all,...
    K, pr_y, N_y, resource_y1, resource_y2, y_medi, x_fc, pr_fc, ...
    rho, crra, sigma_c1, beta_c, beta_f, T1, B1_c, B2_c, B1_f, B2_f, delta, ...
    lr_target, N_f, entry_exitflag_vec, id_ist_old, myopic, Er, cost_weight, use_cost_weight,...
    root_dir)

scale_vec_geq1 = sort(scale_vec(scale_vec>=1),'ascend');
scale_vec_l1   = sort(scale_vec(scale_vec<1),'descend');
scale_vec_custom = [scale_vec_geq1; scale_vec_l1]; % NOT sorted
n_geq1 = length(scale_vec_geq1);

disp(' ')
disp('====================================================')
if medi_cf==0 %
    disp('CF rate adjustment cost under baseline Medicaid')
    disp(' ')
    disp(['y_medi = ', num2str(y_medi)])
elseif medi_cf==1
    disp('CF rate adjustment cost under generous Medicaid')
    disp(' ')
    disp(['y_medi_cf (more generous) = ', num2str(y_medi)])
end

% loop over cost multipliers
for ss_custom = 1:length(scale_vec_custom)

    scale_cost2_1 = scale_vec_custom(ss_custom); % cost multiplier
    ss = find(scale_vec==scale_cost2_1); % position in original scale_vec

    tic0=tic;
    disp(' ')
    disp('**********************************************')
    disp(['ss_custom = ', num2str(ss_custom)])

    % change the cost 
    cost2_1_all_cf = cost2_1_all*scale_cost2_1; % variable cost changes
    cost2_0_cf = cost2_0; % fixed cost held constant

    % initial guess on eqm outcomes
    if ss_custom==1 % baseline estimates
        % use observed data
        p1_ini = D_input.p1; % NN x 1
        p2_ini = D_input.p2; % NN x KM
        n_ini_unmerge = D_input.n_insurers; % NN x 1: 1 for majors, # fringes for the fringe
        n_ini = zeros(N.st,1); % N_st x 1. bring down to market level
        for mm=1:N.st
            n_ini(mm) = unique(n_ini_unmerge(D_input.mkt==mm & D_input.major==0));
        end

    elseif ss_custom==n_geq1+1 % largest value among scales <1
        % use solution from scale=1
        ss_ind = find(scale_vec==1);
        p1_ini_temp = OO(ss_ind).p1_cf; % N.ist x 1 
        p2_ini_temp = OO(ss_ind).p2_cf; % N.ist x KM
        n_ini = OO(ss_ind).n_cf;  % N.st x 1, stores # fringes in each market

        % for prices, change to unmerged version
        p1_ini = unmerge_fringe(p1_ini_temp, D_major, arg_share, D_nF, N_st); % NN x 1
        p2_ini = unmerge_fringe(p2_ini_temp, D_major, arg_share, D_nF, N_st); % NN x KM

    elseif (ss_custom>=2 && ss_custom<= n_geq1)
        % use eqm outcomes from closest scale solved
        p1_ini_temp = OO(ss-1).p1_cf; % N.ist x 1 
        p2_ini_temp = OO(ss-1).p2_cf; % N.ist x KM
        n_ini = OO(ss-1).n_cf;  % N.st x 1, stores # fringes in each market

        % for prices, change to unmerged version
        p1_ini = unmerge_fringe(p1_ini_temp, D_major, arg_share, D_nF, N_st); % NN x 1
        p2_ini = unmerge_fringe(p2_ini_temp, D_major, arg_share, D_nF, N_st); % NN x KM

    elseif ss_custom>=n_geq1+2
        % use eqm outcomes from closest scale solved
        p1_ini_temp = OO(ss+1).p1_cf; % N.ist x 1 
        p2_ini_temp = OO(ss+1).p2_cf; % N.ist x KM
        n_ini = OO(ss+1).n_cf;  % N.st x 1, stores # fringes in each market

        % for prices, change to unmerged version
        p1_ini = unmerge_fringe(p1_ini_temp, D_major, arg_share, D_nF, N_st); % NN x 1
        p2_ini = unmerge_fringe(p2_ini_temp, D_major, arg_share, D_nF, N_st); % NN x KM

    end

    % solve the counterfactual eqm
    [p1_cf, p2_cf, n_cf, ...
        profit_cf, non_converge_all_cf, solve_market_all_cf, fsolve_all_cf, ...
        EV_y_cf, EV_cf, enrollment_cf, enrollment_m_cf, enrollment_f_cf, enrollment_per_m_cf, enrollment_per_f_cf, non_converge_cf, solve_market_cf,...
        s1_mkt_cf, s2_mkt_cf, rs_share_cf,...
        lapse_model_mkt_cf, fsolve_p2_all_cf]...
        = solve_eq_lapse(lapse_u, alpha, ...
        Delta_xi_all, cost1_e_all, cost2_1_all_cf,  mu_f_all, sigma_f_all,...
        cost2_0_cf, D_input, N, ...
        K, pr_y, N_y, resource_y1, resource_y2, y_medi, x_fc, pr_fc, ...
        rho, crra, sigma_c1, beta_c, beta_f, T1, B1_c, B2_c, B1_f, B2_f, delta, ...
        lr_target, N_f, entry_exitflag_vec, id_ist_old, myopic, Er,...
        p1_ini, p2_ini, n_ini, cost_weight, use_cost_weight);

    % save all outcomes
    OO(ss).p1_cf=p1_cf;
    OO(ss).p2_cf=p2_cf;
    OO(ss).n_cf=n_cf;
    OO(ss).profit_cf=profit_cf;
    OO(ss).non_converge_all_cf = non_converge_all_cf;
    OO(ss).solve_market_all_cf = solve_market_all_cf;
    OO(ss).fsolve_all_cf = fsolve_all_cf;
    OO(ss).EV_y_cf = EV_y_cf;
    OO(ss).EV_cf = EV_cf;
    OO(ss).enrollment_cf = enrollment_cf;
    OO(ss).enrollment_m_cf = enrollment_m_cf;
    OO(ss).enrollment_f_cf = enrollment_f_cf;
    OO(ss).enrollment_per_m_cf = enrollment_per_m_cf;
    OO(ss).enrollment_per_f_cf = enrollment_per_f_cf;
    OO(ss).non_converge_cf=non_converge_cf;
    OO(ss).solve_market_cf=solve_market_cf;
    OO(ss).s1_mkt_cf = s1_mkt_cf;
    OO(ss).s2_mkt_cf = s2_mkt_cf;
    OO(ss).rs_share_cf = rs_share_cf;
    OO(ss).lapse_model_mkt_cf = lapse_model_mkt_cf;
    OO(ss).fsolve_p2_all_cf = fsolve_p2_all_cf;

    toc0 = toc(tic0);
    disp(' ')
    disp(['Time elapsed for current ss_custom (min) = ', num2str(toc0/(60), 5)]);
end

% solve_market=1 when cost estimates make sense 
ssolve_market_all = OO(1).solve_market_all_cf; % N.ist x 1
ssolve_market     = OO(1).solve_market_cf; % N.st x 1

% check convergence 
nnon_converge_all = 0; % N.ist x 1
nnon_converge = 0; % N.st x 1
for ss=1:length(scale_vec)
    nnon_converge_all = nnon_converge_all + OO(ss).non_converge_all_cf;
    nnon_converge = nnon_converge + OO(ss).non_converge_cf;
end

if medi_cf==0
    ss_ind = find(scale_vec==1);

    p1_eq            = OO(ss_ind).p1_cf; % N.ist x 1
    p2_eq            = OO(ss_ind).p2_cf; % N.ist x KM
    n_eq             = OO(ss_ind).n_cf;  %N.st x 1
    profit_eq        = OO(ss_ind).profit_cf; % N.ist x 1
    non_converge_all = OO(ss_ind).non_converge_all_cf; % N.ist x 1
    solve_market_all = OO(ss_ind).solve_market_all_cf; % N.ist x 1
    fsolve_all       = OO(ss_ind).fsolve_all_cf; % N.ist x 1
    EV_y             = OO(ss_ind).EV_y_cf; % N.st x N_y
    enrollment_eq    = OO(ss_ind).enrollment_cf; % N.st x 1
    enrollment_m_eq  = OO(ss_ind).enrollment_m_cf; % N.st x 1
    enrollment_f_eq  = OO(ss_ind).enrollment_f_cf; % N.st x 1
    non_converge     = OO(ss_ind).non_converge_cf; % N.st x 1
    solve_market     = OO(ss_ind).solve_market_cf; % N.st x 1
    rs_share_eq      = OO(ss_ind).rs_share_cf;  % N.ist x 1

    [fit_off_st, fit_off_ist] = report_baseline(enrollment_eq, non_converge, solve_market,...
        D, K, N, D_enroll_st);

elseif medi_cf==1
    fit_off_st  = zeros(N.st,1);
    fit_off_ist = zeros(N.ist,1);
end

% compute changes
disp(' ')
disp('====================================================')
if medi_cf==0
    disp('Baseline Medicaid')
elseif medi_cf==1
    disp('Generous Medicaid')
end

ss_ind = find(scale_vec==1);
ns = length(scale_vec);
RS = report_outcome(OO, ns, ss_ind, ...
    ssolve_market_all, nnon_converge_all, ssolve_market, nnon_converge,...
    D, K, N, D_enroll_st, pr_y, alpha, B1_c, beta_c, T1, B2_c,...
    fit_off_st, fit_off_ist);

% save
if medi_cf==0
    save('temp_variable_cost.mat')
elseif medi_cf==1
    save('temp_variable_cost_medi_cf.mat')
end

% Figure A.3
folder_tex_cf = '...';
xl = '$c^1_{jk}$ multiplier'; % x-axis label
fig_title = 'Counterfactual rate stability regulation - endogenous lapses';
lr = 0;
export_result_lapse(RS, scale_vec, xl, fig_title, lr, lr_target, folder_tex_cf);
close all

clear OO

